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Summary 

This paper describes a method to solve a system of 
N linear equations in N steps. A quadratic form is devel- 
oped involving the sum of the squares of the residuals of 
the equations. Equating the quadratic fonn to a constant 
yields a surface which is an ellipsoid. 

For different constants, a family of similar ellipsoids can 
be generated. Starting at an arbitrary point an orthogonal 
basis is constructed and the center of the family of similar 
ellipsoids is found in this basis by a sequence of projec- 
tions. The coordinates of the center in this basis are the 
solution of the linear system of equations. A quadratic 
form in N variables requires N projections. That is, the 
current method is an exact method. It is shown that the 
sequence of projections is equivalent to a special case of 
the Gram-Schmidt orthogonalization process. The current 
method enjoys an advantage not shared by the classic 
Method of Conjugate Gradients. The current method can 
be extended to nonlinear systems without modification. 
For nonlinear equations the Method of Conjugate 
Gradients has to be augmented with a line-search 
procedure. Results for linear and nonlinear problems 
are presented. 

Introduction 

The computation method described herein belongs to the 
general class of methods where the iterates are derived 
from information available about the function and its 
derivatives at each step. It applies to functions where the 
minimum value of the function is known to be zero. In 
order to apply the method to the general multidimensional 
root finding problem, a function is introduced which is 
the sum of the squares of the individual functions. This 
function is called the squared-error function. It is positive 
definite and has a global minimum of zero at all solutions 
of the original set of equations. 

Now consider the case of linear equations. In this case the 
sum of the squares of the residuals is considered. The 
squared-error function is a quadratic function of the 
N variables. This quadratic form is positive definite and 
has a global minimum of zero at all solutions of the 
original set of equations, provided that the rank of the 


coefficient matrix is N. Equating the quadratic form to a 
constant yields a surface which is an ellipsoid. A basis for 
this space is formed from the set of vectors reciprocal to 
the conjugate radii of the ellipsoid. The conjugate radii and 
the reciprocal vectors are generated by an interlocking 
iterative procedure. This procedure is followed to solve a 
system of linear equations. The solution of N linear 
equations is accomplished in N iterations. That is, the 
current method is an exact method for linear systems. 

An advantage of this method is that it can be extended to 
systems of nonlinear equations without modification. The 
Method of Conjugate Gradients (ref. 1) solves a linear 
system of dimension N in N steps. However, in order to 
extend this method to nonlinear equations a line search 
procedure has to augment the original method. The line 
search procedure is required at each stage to generate 
successive iteration points. At each stage the step size is 
the solution of a minimization problem along a line in a 
selected direction (ref. 2). 

This paper has a twofold purpose. One purpose is to 
present the method of conjugate radii for linear systems. 
Another purpose is to extend the method to nonlinear 
systems. “Nonlinear problems can in most cases be solved 
only by approximating them by linear problems” (ref. 3). 
The approach of the current method to solve a nonlinear 
system is to approximate it by a quadratic function. 

The current method is quadraiically convergent That is, 
for a quadratic function of N variables, the root will be 
found in N iterations. In the nonlinear case, assume the 
squared-error function can be expanded in a Taylor series 
about the root. The matrix whose elements are the second- 
order partial derivatives of the function is called the 
Hessian matrix. It is symmetric and positive definite. If it 
is assumed that the function can be approximated near the 
root by a quadratic form, the elements of the Hessian 
matrix are the coefficients of the quadratic form. For 
nonlinear functions, as the iterates approach the root, the 
squared-error function is more closely quadratic and hence 
convergence is more nearly quadratic. 

After the relevant properties of central quadratic surfaces 
are identified, the method is described for linear systems 
and an algorithm is formulated. A linear and a nonlinear 
example are treated by the same algorithm. The relation 


between the Gram-Schmidt orthogonalization procedure 
and the current method is discussed, as is the duality 
between the current method and the Method of Conjugate 
Gradients. 

Method of Conjugate Radii 

Problem Statement 

The system of N linear equations A • X = b is solved by 
introducing the function 

2/ = (A-x-b)-(A-x-b) = x-A’ A x-2x A' b + b b 

0 ) 

The raised dot in an expression such as A • X denotes 
matrix multiplication. Square matrices appear as bold 
capital letters and vectors appear as bold lower case letters. 
In the dot product, vectors that appear before the dot are 
understood to be row vectors and those that appear after are 
understood to be column vectors. The transpose of a 
matrix is denoted by a prime. 

For an arbitrary x the function (1) is the square of the 
magnitude of the residual b — A * X which is a measure 
of the error incurred by the choice of x. In carrying out 
the computations, it is necessary to evaluate the gradient. 

The formula for the gradient is 

g = V/ = A'-(Ax-b) (2) 

When the squared error is held constant, equation (1) 
represents a family of quadratic surfaces. Each constant 
value leads to a different surface. For a system of two 
linear equations the significance of the quadratic surfaces 
can be visualized. If the equations are expressed in terms 
of the usual Cartesian coordinates (x,y), and if the squared 
error represents altitude above the (x,y) plane, a three- 
dimensional surface is obtained. The surface has the form 
of a bowl, an elliptic paraboloid with vertex at the point c 
in the (x,y) plane. Figure 1 shows multiple level lines of 
such a surface. Two of these level lines are projected onto 
the (x,y) plane. The vector c is the solution of the system 
of equations, that is, A • C = b . The projections of the 
level lines of the squared-error function onto the (x,y) 
plane represents a family of concentric ellipses, quadratic 
surfaces, similar to each other. These are called error- 
ellipses (ref. 3). In higher dimensions, an error-ellipsoid is 
obtained. 

The problem of solving the linear system A • X = b for 
x is equivalent to the problem of determining the center c 
of the family of similar error-ellipsoids generated by the 
squared-error function. 


In view of the equivalence of the two problems, the 
following problem is addressed herein. The problem to be 
solved is: find the origin of a system of coordinates in 
which a family of similar central quadratic surfaces is 
embedded starting at an arbitrary point. 

The central quadratic surfaces associated with the problem 
A • r — 0 are related to the non-central quadratic surfaces 
associated with the problem A * X = b by means of the 
coordinate transformation r — X — C . It is to be empha- 
sized in what follows that x and r are the same point, and 
this point has different labels and components, in the two 
coordinate frames. 

The point x = c corresponds to the center of the family of 
similar error-ellipsoids in the x coordinate frame, and the 
point r = 0 corresponds to the center of the same family 
in the r coordinate frame. It will be shown that the 
sequence of operations that leads to the origin in the r 
coordinate frame, starting from an arbitrary point r, leads 
to the solution c in the x coordinate frame, starting from 
an arbitrary point x. 

The properties of central quadratic surfaces play a 
significant role in the analysis which will be presented. 
Denoting A - A by H, the expressions for the squared- 
error function and the gradient can be written in the r 
coordinate frame as follows: 


2/ = r H r 

(3) 

g = Hr 

(4) 


The representation of planes in space and the essential 
properties of central quadratic surfaces which will be used 
in the analysis are identified. The three-dimensional case 
will be discussed. 

Representation of Planes and Properties of 
Central Quadratic Surfaces 

In analytic geometry the intercept form of a plane is 
written as 

r n = 1 

where the components of n are the reciprocals of the 
intercepts of the plane on the coordinate axes and r is the 
position vector of a general point on the plane. 

This equation when rewritten as follows reveals some 
interesting features: 

|r|cos(n,r)|n| = 1 

Here |r| and |n| denote the magnitudes of the vectors r 
and n respectively and (n,r) denotes the angle between n 
and r. Now |n| , the length of n, is the reciprocal of the 
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perpendicular distance from the origin to the plane since 
|r|cos(n,r) is that distance. As stated, the components 
of n are the reciprocals of the intercepts of the plane on 
the axes. 

Gibbs (ref. 4) points out that the vector n can be used to 
denote the plane. It may be taken as the vector coordinate 
of the plane. To designate a plane denoted by the vector n 
the notation (n,l) is used. This is the notation of analytic 
geometry to specify the “coordinates of a plane.” 

If the elements of H are constant, then 

r • H • r = const 


is a quadratic in r . It represents an ellipsoid with center 
at the origin. 


Hence for the problem at hand, the following family of 
similar ellipsoids is considered: 


r H r 

2 / 


= 1 


(5) 


Define the normalized gradient 




( 6 ) 


reciprocal of the perpendicular distance from the origin to 
the tangent plane at . The plane coordinates of the 
tangent plane are ( n v 1). 

n i 

The vector is the vector drawn from the center, 

n i n i 

the origin of the r coordinate frame, perpendicular to the 

tangent plane (ref. 4). See figure 2. The standard approach 
for solving a system of equations is to assume an initial 
approximation T x and to proceed to an improved approxi- 
mation r 2 by using an iterative formula of the form 


r 2 = r i”‘ 


n, n. 


The significance of the vector drawn from the origin 
perpendicular to the tangent plane is that it represents, in 
the terminology of the standard approach, the negative of 
the “direction of search,” and its magnitude determines the 
“step size” in this direction. As explained below, the 
strategy of establishing a plane at a point I* and moving 

1 

toward the origin a distance - — - is the basis of the 

hi 

current method. 


From this definition it follows that n • r = 1 . 

Figure 2 illustrates the foregoing relations in two 
dimensions. The tangent plane to an ellipse at the point R 
is shown. The normalized gradient at this point is also 
shown. This is the vector originating at the point R and 
terminating at the point N. The center of the ellipse is 
labeled C and this point is the origin of the r coordinate 
frame. The point labeled 0 in this figure is the origin of 
the x coordinate frame. The components of the vector n, 
which is the normalized gradient referred to the origin in 
the r coordinate frame, are (1,1). The intercepts of the 
tangent plane when referred to the same origin are (1,1). 
The vector n represents the tangent plane. The length of 
this vector is V2 . The distance of the tangent plane 

from the origin is l/ V2 which is the reciprocal of the 
length of n. Also, the vector which represents the plane 
tangent to the similar inner ellipse at the point (0,.5) with 
intercepts (.5, .5) would have components (2,2) and length 
2 ■ V2 . The distance of this plane from the origin is 
l/ (2 • a/ 2) . That is, the nearer a plane is to the origin 
the longer the vector will be that represents it. 

At any point r } of the quadratic surface (5), n p the vector 
given by equation (6) at that point, is normal to the 
tangent plane at . The equation of the tangent plane at 
Tj is n a • (r — r } ) = 0. Since Tj • llj = 1, this can be 
written n x • r = 1. Now the length of llj is the 


Any two vectors r and S are referred to as conjugate 
vectors if r H s =0 . The tangent plane at the point 
Tj is (Iij, l).The plane through the origin parallel to 
this tangent plane has plane coordinates ( llj , 0). It is 
referred to as the diametral plane conjugate with T x . 

Any vector r in that plane satisfies the relation 

r H Tj = r nj =0. 

If three points located on similar ellipsoids of the form (5) 
satisfy the equations 

r, • H • r 2 = r 2 • H • r 3 = r 3 • H - r, =0 

the position vectors (ii,r 2f r 3 ) are said to form a 
conjugate set. The vectors (r,,r 2 ,r 3 ) and (npn 2 ,n 3 ) 
form reciprocal sets. See reference 4. The relation 
ij H r 2 is symmetric in the indices 1 and 2; that is, 

Tj * H • r 2 = r 2 • H • T x . Hence, if llj * r 2 = 0 then 

n 2 • r i = 0 . ‘ 

Method of Solution 

The method of solution is similar to the process of 
assigning Cartesian coordinates (x,y,z), given the 
Cartesian frame, to a point P in space. The frame is given 
by giving three mutually perpendicular coordinate planes. 
Their lines of intersection are the axes of the coordinates 
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and are called the x axis, y axis, and z axis. The point of 
intersection of the coordinate planes is the origin. Let P 
be any point in space and let three planes be drawn 
through P perpendicular to the coordinate planes and 
cutting the axes at x, y, and z. These three planes 
together with the coordinate planes bound a rectangular 
parallelepiped, of which P and the origin 0 are opposite 
vertices. The three edges x, y, and z are called the 
rectangular coordinates of P. That is, the rectangular 
coordinates of P are equal to its perpendicular distances 
from the coordinate planes. 

The method of solution described herein follows the 
procedure described in the previous paragraph starting at 
the point P. That is, the origin is located with respect to 
P. The length of the edges of the parallelepiped are 
conveniently found using the vector representation of a 
plane. The strategy at each iteration point for the construc- 
tion of a plane is enabled by the relations satisfied by the 
mutually conjugate position vectors generated up to that 
point. 

The construction will be carried out for the three- 
dimensional case. Figure 3 shows the parallelepiped. The 
vertex of the parallelepiped opposite the origin 0 is labeled 
R x . The vertex of the parallelogram in the plane of the 
three points 0, R 2 , and R 3 opposite the origin 0 is 
labeled . The vertex of die line connecting the points 0 
and R 3 opposite the origin 0 is labeled R 3 . 

In three dimensions, three iteration points (r p r 2 ,r 3 ) will 

be generated. It will be demonstrated that the iteration 
points form a conjugate set. In carrying out the procedure 
the following results will be demonstrated: 

1) The vectors n p llj + n 2 ,nj + n 2 + n 3 form a 
three-dimensional orthogonal basis. 

2) The vectors reciprocal to the vectors 

Ilpllj + n 2 ,n 1 + n 2 + n 3 satisfy the following 

relations 


_ r + n 2 

3 (n 1 +n 2 )-(n, +n 2 ) 

°i + n 2 + n 3 

(iij + n 2 + n 3 ) • (nj + n 2 + n 3 ) 


(7) 

( 8 ) 
(9) 


The reciprocal vectors are the right hand members of 
equations (7), (8), and (9) respectively. 

Carrying out the steps of the procedure will be a 
constructive demonstration of the above results. The 


procedure can be regarded as a coordinate transformation 
from a given basis to a new basis. The coordinate 
transformation will be accomplished by a sequence 
of projections orthogonal to each of the vectors 
llpllj + n 2 + n 3 . These vectors themselves 

are obtained during the sequence. 

The problem of finding the origin will be carried out in 
the new basis. Referring to figure 3, R { will be projected 
onto the parallelogram yielding . Then R 1 will be 
projected onto the line 0 - R 3 yielding R 3 . Finally, R 3 

will be projected onto the origin. These projections will 
be identified at each step. 

At each step of the iteration, it is necessary to compute 
the components of the original basis n p n 2 ,n 3 in terms 

of the new basis n^llj + n 2 ,nj + n 2 + n 3 . These 

components are needed to carry out the next step of the 
iteration or to show that the iteration is complete. 


The process starts at an arbitrary point. In general, once r 
is known, n can be calculated using (6). At the point r } , 
construct the plane (llj , 1). That is, use the point 
normal form of the equation of a plane to write 
llj • (r — r t ) = 0. Since llj • r 2 = 1, the intercept form 
of the plane is I1 1 • r = 1 . This plane is the tangent plane 
at the point r x . Then construct a plane parallel to the 
plane at that passes through the origin. This is the 
plane (lip 0). This is the first coordinate plane. The 
distance of the original plane (n l , 1) from the origin is 


n, 


■ . The location of a point r 2 that lies on the 


Cartesian coordinate plane ( lip 0), the parallelogram in 
figure 3, is obtained by projecting onto that plane. 
That is 


However, lies on the tangent plane (llj ,1); therefore, 
llj • I\ = 1. Hence, the location of the point is given by 
the equation 

n i 

r 2 = i\ i— (10) 

n i n i 

The scalar product I1 1 ■ r 2 = • H • r 2 = 0 shows that 

the two position vectors I\, r 2 form a conjugate set since 
n i * r 2 “ 0 - By symmetry n 2 * T y = 0 . That is, 
n i r 2 = Tj H*r 2 = r, n 2 . 
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Hence, the vectors (r p r 2 ) and (ll],n 2 ) form reciprocal 
sets. 

Therefore 

(n 1 +n 2 )(r 1 -r 2 ) = 0 (ll) 

Equations (7) and (11) show that 

n, (n t + n 2 ) = 0 (12) 

The following relation follows from equation (12): 

n 2 -(n! +n 2 ) = (n 1 +n 2 )-(n 1 +n 2 ) (13) 

The relations (12) and (13) are needed to express the 
components of the basis n 1 ,n 2 ,n 3 in terms of the new 
basis Hpllj +n 2 +n 3 . These relations for 

two dimensions are summarized in the following 
multiplication table: 



F 1 _I 2 


* 

D I 

n,+n 2 

2 

**i 

2 

n, +n 2 

n l 

1 

0 

n 2 

-1 

1 


In the table an entry in the dashed row is equal to the entry 
in the corresponding column in the row directly below it. 

The projection of n 2 onto the (n p 0) plane is llj + n 2 . 
That is 

11 , 11 ] 

n, + n 2 = n 2 ^ • n 2 (14) 

H] n. 

The multiplication table was used to arrive at this 
equation. 

Equation (12) suggests the vector (n, + n 2 ) can be used 

as the normal to construct a plane perpendicular to the 
plane ( n p 0). At the point r 2 , construct the plane 
0*i +n 2 ,l). The plane parallel to this plane that passes 
through the origin is (n, + n 2 ,0). Since (n p 0) is 
the first coordinate plane, it is appropriate to select 
(n, + n 2 ,0) as the second coordinate plane. 

The location of a point r 3 that lies on the second 
coordinate plane, (II] + n 2 , 0), is obtained by projecting 
r 2 onto that plane. That is 

_ _ (n, +n 2 )(n 1 +n 2 ) _ 

r, = r 2 r 2 


However, r 2 lies on the plane (ll, + n 2 , 1) ; therefore, 
(n, +n 2 )-r 2 = 1. 


Hence, the location of the point is given by the equation 


r 3 = r 2 - 


n i +11 2 


(n,+n 2 )-(n,+n 2 ) 


(15) 


Using equation (15) to represent r 3 and utilizing the 

relations (12) and (13) which are summarized in the 
multiplication table, it is readily verified that 
It] * r 3 = n 2 * r 3 = 0 and by symmetry 

n 3 -r, =n 3 r 2 =0. 


The geometric interpretation of this result is interesting. 
The equation n, * r = 0 is the equation of the diametral 
plane conjugate with Tj , and the equation n 2 • r = 0 is 
the equation of the diametral plane conjugate with r 2 . 
Simultaneously, these two equations represent the line of 
intersection of the two diametral planes. The results, 
n, * r 3 = n 2 ■ r 3 = 0 , establish that r 3 is along this line 

of intersection. This is the line 0-R 3 shown in figure 3. 

The conjugate relation is symmetric. Hence, r, is 
in the plane conjugate with r 3 as is r 2 . That is, 
r 2 • H * r 3 = r 3 * H * Tj = 0 . This result along with 
the result r, * H * r 2 =0 established previously indicates 
that the position vectors (r, , r 2 , r 3 ) form a conjugate set. 

Since the vectors (r],r 2 ,r 3 ) and (n],n 2 ,n 3 ) form 


reciprocal sets, it follows 

(n, + n 2 + n 3 ) • (r 2 - r, ) = 0 (16) 

(**1 **3) ‘ (**3 — **2 ) = ® ( 17 ) 

Equations (7) and (16) show 

n, -(n, +n 2 + n 3 ) = 0 (18) 

Equations (8) and (17) show 

(it] +n 2 )-(n 1 + n 2 +n 3 ) = 0 (19) 

From equations (18) and (19) it follows 

n 2 (n, + n 2 + n 3 ) = 0 (20) 


The final relation needed to express the components 
of the basis Il],n 2 ,n 3 in terms of the new basis 
ll,, 11] + n 2 ,n] + Il 2 + n 3 can be obtained by combining 
equations (18) and (20). This relation is 

n 3 -(ii] +n 2 +n 3 ) = (n, +n 2 +n 3 )*(nj +n 2 +n 3 ) 

( 21 ) 
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The foregoing relations for three dimensions are used to 
update the two-dimensional multiplication table. The 
current relations are summarized in the following 
multiplication table: 



r i -I 2 

r 2 -r 3 


* 


n, +n 2 

iij +n 2 +n 3 

2 

n, 

2 

iij +n 2 

2 

n 3 + n 2 + n 3 

n . 

1 

0 

0 

n 2 

-1 

1 

0 

n 3 

0 

-1 

1 


The projection of n 3 onto the plane ( Iij + n 2 , 0) is 
n 1 + n 2 + n 3 . That is 


n, + n 2 + n 3 = n 3 


(n t +n 2 )(n, +n 2 ) p 
(n, +n 2 )-(n, +n 2 ) 

( 22 ) 


The current multiplication table was used to arrive at this 
equation. 

Equation (19) suggests the vector (ilj + n 2 + n 3 ) can be 

used as the normal to construct a plane perpendicular to 
the plane ( n 3 4- n 2 , 0). At the point r 3 , construct the 

plane (iij + n 2 4- n 3 ,1). The plane parallel to this plane 
that passes through the origin is (ti x 4- n 2 4- n 3 ,0). 

Equation (18) shows that this plane is perpendicular 
to the plane (n 1? 0), and equation (19) shows that it is 
perpendicular to the plane ( Ilj 4- n 2 , 0). Since (n l -f n 2 , 
0) is the second coordinate plane, it is appropriate to select 
(ilj 4- n 2 4- n 3 ,0) as the third and final coordinate plane. 


The location of a point r 4 that lies on the third coordinate 
plane (Ilj + n 2 4- n 3 , 0) is obtained by projecting r 3 
onto that plane. That is 

r ( n ! + n 2 + n 3 )(n 1 + n 2 + n 3 ) f 

4 3 (iij +n 2 4-n 3 )-(n r + n 2 +n 3 ) 3 

However, r 3 lies on the plane (ilj 4- n 2 4- Il 3 ,1) ; there- 
fore, (Ilj 4- n 2 4“ n 3 ) • r 3 = 1 . Hence, the location of the 
point is given by the equation: 

n, +n 7 +n. 

r — r 1 2 (23) 

4 3 (iij 4-n 2 +n 3 )*(n 1 +n 2 + n 3 ) 


This equation will be used to show r 4 = 0 . 

In order to express a vector in a basis (v x , r 2 , r 3 ) , its 
components in the reciprocal basis (n 1 ,n 2 ,Il 3 ) are 
formed. Hence r 4 can be expressed as follows: 

r 4 = (r 4 • n , )r, + (r 4 • n 2 )r 2 + (r 4 • n 3 )r 3 (24) 

Using equation (23) to represent r 4 , the components can 
be calculated using equations (17), (19), and (21) which 
are summarized in the current multiplication table. All the 
components vanish, hence r 4 = 0 . 

The above procedure, carried out in the r coordinate frame, 
led to the center, r = 0, starting at an arbitrary point. The 
same procedure, carried out in the x coordinate frame, will 
also lead to the center, x = c, starting at an arbitrary 
point. That is, r 4 = 0 corresponds to x 4 = C . 


Algorithm 

For any two points, 1 and 2 for instance, it is the case 
that Xj — X 2 — T x — r 2 . This equation is correct since it 
is independent of the origin. Using this fact and express- 
ing the terms involving r in equations (7), (8), and (9) in 
terms of x and rearranging those equations leads to the 
following procedure in the x coordinate frame. 

The procedure in the x coordinate frame starting from an 
arbitrary point, X v can be summarized as follows: 


X 3 = X 2 - 


n x +n 2 


c = x 3 — 


(n,+n 2 )-(n,+n 2 ) 

n,+n 2 +n 3 

(n, + n 2 + n 3 ) • (n, + n 2 + n 3 ) 


Given Xj the above equations can be evaluated in 
sequence using the relation n = g/2/. In general, 
once x is known, g and 2f can be calculated using 
equations (1) and (2) for the linear case. For the linear 
case, it is understood that A is an NxN matrix of full 
rank. For the nonlinear case the gradient vector g = V/ 
is obtained by calculating the required partial derivatives 
of the squared-error- function 2f. 

For N dimensions the algorithm is implemented by 
carrying out the following steps: 
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Set n = 0 

Set x = 0 or any arbitrary value 

for index = 1 : 1 :N where N is the number of 
unknowns 

n = n + g/2f 

n 

X = X 

n n 

end 

return x 

Examples 
A Linear Example 

A two-dimensional example will be presented. This 
example will illustrate the application of the method and 
will also illustrate the relation between the basis (r 19 r 2 ) 
and its reciprocal (llj ,n 2 ) and the relation between the 
orthogonal basis (llj ,nj + n 2 ) and its reciprocal 
, n x n x +n 2 

( , 7 r). The solution of the 

n i n 2 (nj +n 2 )-(n 1 +n 2 ) 
following linear system is considered: 

3x+y = 5 
Jt+2 y - 5 

Equations (1 ) and (2) are used to compute 2f and g . 

Figure 4 shows the projection of two level lines of the 
squared-error function, 2f, onto the (x,y) plane. The 
solution of the system of equations is (x,y) = (1,2). 

This point is labeled C in the figure. The origin of the 
x coordinate frame is labeled 0. The starting point was 
(x,y) = (1,3). The two vectors forming the basis (r l9 r 2 ) 
are labeled (Rj,R 2 ) in the figure. The vectors (n,,n 2 ) 
reciprocal to them are also shown and are labeled (Nj,N 2 ). 
In addition, the parallelogram of which Rj is a diagonal is 

shown. The vertices of this parallelogram are labeled C, 

P, Rj, and R 2 . The portion of the parallelogram that is not 

obscured by solid lines is shown as the dashed line. The 
line Rj-P is tangent to the ellipse at the point Rj. The 

vector Ilj + n 2 terminates at the point Q. The magnitude 
of both of the vectors and I1 1 + Il 2 is V2. 

✓ Hi "t" n 2 

The vectors ( , -) reciprocal 

nj n 2 (n, +n 2 )-(n, +n 2 ) 

to the base vectors (ilj ,Itj + n 2 ) are not labeled in the 

figure. They are in the same direction as the base vectors 


and their magnitudes are the reciprocals of the base 
vectors, namely l/ V2 . The first reciprocal vector 
terminates at the point P and the second terminates at the 
point R 2 . They are the edges of the parallelogram of 

which Rj is the diagonal. 

In order to describe the sequence of projections, the same 
problem as above is discussed but with a different starting 
point. Figure 5 identifies the projection planes. All the 
labels shown in this figure have the same significance as 
the labels in figure 4. The starting point used in figure 5 
was (x,y) - (2,0). The identification of the four planes 
used in the construction is as follows: The (n l9 l) plane 
is the horizontal plane that passes through the point . 
The parallel horizontal plane (n p 0) is the plane that 
passes through the point C, the center. The (lt 1 +n 2 ,l) 
plane is the vertical plane that passes through the 
point /?2 . The parallel vertical plane (lij + n 2 ,0) is the 
plane that passes through the point C. 

The projections described earlier were performed. First, 

T x which is located in horizontal plane (n 19 l) was 
projected onto the parallel (lij ,0) plane, as was n 2 . 

This led to the position vector r 2 and the vector 
Ilj + n 2 respectively. The vectors r 2 and Itj + Il 2 
both terminate at the point . Second, r 2 , which is 
located on the vertical plane (ilj + n 2 , 1) , was 
projected onto the parallel (ll 1 +n 2 ,0) plane. This 
results in the point C which is the solution of the system 
of equations. The orthogonal planes which were 
constructed are (ilj , 0) and (n x + n 2 , 0) . The vectors 
, n 2 n 1 + n 2 

( , 7 -) which terminate at 

n, n 2 (iij +n 2 )-(n 1 +n 2 ) 

the points (P , ) are the edges of the parallelogram of 

which Rj is the diagonal. 

The results shown in figure 5 permit a comparison to be 
made between the current method and the Method of 
Conjugate Gradients. If the calculation was performed by 
that method starting at the point R x , the next point would 
be the point labeled CG in the figure. This point is 
located by minimizing the squared-error function along the 
line R l ~R 2 shown in the figure. The radii C-CG and 
C-P are conjugate radii. 

A Nonlinear Example 

The extension of the current method to nonlinear systems 
will be illustrated for Rosenbrock’s function (ref. 5). 

1f{x,y) = 100 (y-x 2 ) 2 +(1-*) 2 
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The current method can be applied to find the root 
(jt, y) = (1,1), where f(x,y) = 0, since Rosenbrock’s 
function represents the sum of the squares of two 
individual functions. The gradient vector g = V/ is 
readily found by calculating the required partial derivatives 
of the function f. A contour map is shown in figure 6 
which shows the solution path to the root for the current 
method starting at the point (-1.2,1). The function value 
was reduced to 4 x 10“ 13 in 13 iterations. This compares 
with Fletcher and Reeves’ value of 1 x 10 -8 after 
27 iterations using a line search procedure together 
with the Method of Conjugate Gradients (ref. 2). 

Another contour map is shown in figure 7 which shows 
the solution path to the root for the current method 
starting at the point (-1.92,2). Results for this starting 
point are given in reference 6. These results are compared 
with the current method in the table below. For these 
results the function value was less than 1 x 10 - " 4 . 


Method 

Iterations for 
convergence 

Gauss-Newton 

48 

Levenberg-Marquardt 

90 

Current 

10 


Discussion of Results 


The sequence of projections orthogonal to each of the 
vectors n l9 n x 4* n 29 n { + Ii 2 + n 3 is identical to the 
Gram-Schmidt orthogonalization process applied to the 
vectors n^n^I^. The results of the Gram-Schmidt 
process are shown below: 

Hi = Dj 


111 ~ ^2 — R 2 


n i n i 

"i n, 


n, 


(25) 


n.n, 

n, + n 2 + n 3 = n 3 • n 3 

n r n i 


(n 1 +n 2 )(n 1 +n 2 ) p 
(n, +n 2 )-(n, +n 2 ) 


(26) 


Equations (12) and (18) combined show that • n 3 = 0. 

The same result is obtained by consulting the current 
multiplication table. It can be seen that there is no need to 
perform the first subtraction in equation (26) since the 
term to be subtracted is zero. Hence, the sequence of 
projections is identical to the Gram-Schmidt 
orthogonalization process applied to the vectors 
11] ,n 2 ,n 3 . 


That they are identical for all cases can be seen by 
rewriting equations (25) and (26) employing equations (7) 
and (8) which give alternate expressions for the vectors 
reciprocal to Ilj and lt 1 + n 2 . 

The Gram-Schmidt equations (25) and (26), rewritten 
employing equations (7) and (8), are shown below: 

it, + n 2 = n 2 - n, (r 3 - r 2 ) • n 2 (27) 

n,+n 2 +n 3 =n 3 - n^r, - r 2 ) • n 3 - (n, + n 2 )(r 2 - r 3 ) • n 3 

(28) 

The scalar products involving the vectors of the reciprocal 
sets (ii,r 2 ,r 3 ) and (n 1 ,n 2 ,n 3 ) which are needed to 

evaluate the terms in equations (27) and (28) are arranged 
in a table below: 


(i\ -r 2 )-n 2 


(n -r 2 )-n 3 

(r 2 -r 3 )-n 3 


The only non- vanishing scalar products lie along the 
diagonal. For a given row, the terms that appear before the 
term on the diagonal all involve vanishing scalar products 
of the reciprocal sets of vectors. Hence, only one subtrac- 
tion per step is required to carry out the Gram-Schmidt 
process, and the term subtracted at each step is identical to 
the term subtracted in the projection process. The same 
result is obtained directly using equations (25) and (26) if 
the current multiplication table is consulted to determine 
the value of the scalar products. 

It should be noted that the Gram-Schmidt procedure can be 
used to derive the basic algorithm of the Method of 
Conjugate Gradients (ref. 7). As is the case here, in that 
derivation only one subtraction per step is required. 

There is a relatively direct duality between the Method of 
Conjugate Gradients and the current method. In all 
approaches to root finding, starting at an arbitrary point, 
the following requirements have to be addressed: In 
proceeding from one iteration point to the next, a step size 
and a direction are required. In the Method of Conjugate 
Gradients, consecutive directions are in conjugate 
directions. The step size is such that the gradients at 
consecutive points are orthogonal. In the current method, 
consecutive directions are orthogonal. The step size is 
such that consecutive points lie on conjugate radii. 
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For the three-dimensional case, it is of interest to note 
that the calculation performed herein parallels a similar 
calculation performed in x-ray diffraction studies. The 
vectors in the r coordinate frame correspond to lattice 
vectors, and the vectors in the n coordinate frame corre- 
spond to reciprocal lattice vectors. In x-ray diffraction 
studies, the purpose of the calculation is to determine the 
spacing between lattice planes. This interplanar spacing is 
used to determine the condition for scattering-in-phase 
known as Bragg’s law. 

For systems of nonlinear equations, the current method is 
a general method that converges quadratically to the 
solution. That is, for a quadratic function of N variables, 
it converges to the solution in N steps. The current 
method enables a nonlinear system to be approximated 
by a quadratic function. 
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Figure 1. The squared-error function. 
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